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Coronal  Heating:  by  Means  of 
Helical  Magfnetohydrodynamic  Turbulence 

1.  Introduction 

This  paper  reports  the  results  of  numerical  simulation  of  the  viscoresistive  decay 
of  helically  turbulent  magnetofluids.  This  decay  process  has  recently  come  to  be  of 
interest  because  of  the  possibility  that  it  might  contribute  to  the  heating  of  the  solar 
corona.  It  is  by  now  generally  accepted  that  the  corona  is  heated  by  the  dissipation 
of  magnetic  energy  (see,  c.  g.,  Chapter  6  of  Priest  1982  and  references  therein). 
The  coronal  magnetic  field  is  believed  to  be  continually  perturbed  as  the  footpoints 
of  the  coronal  flux  tubes  aire  convected  about  by  the  turbulent  motions  of  the  nearly 
ideal,  high  0,  photospheric  plasma.  Both  random  displacements  and  helical,  twist¬ 
ing  motions  (perhaps  by  Coriolis  forces  due  to  solar  rotation)  are  possible.  These 
magnetic  field  disturbances  propagate  into  the  corona,  driving  the  plzisma  there 
into  an  excited,  nonequilibrium  state.  As  the  plasma  settles  into  a  new,  quiescent 
state,  the  magnetic  energy  is  released  directly  as  heat  energy  by  Ohmic  dissipation 
and  indirectly  by  viscous  dissipation  of  the  induced  fluid  motions.  The  character¬ 
istic  structure  associated  with  this  dissipation  is  the  familijir  electric  current  sheet 
located  at  a  region  where  the  magnetic  field  changes  sign  {e.g.,  Syrovatskii  197S. 
Matthaeus  and  Montgomery  19S1,  R.  Dahlburg  et  al.  19S6a,  Spicer  et  al.  19S6). 

Only  a  small  portion  of  the  magnetic  energy  in  the  corona  ends  up  being  trans¬ 
formed  into  heat.  The  overwhelming  portion  exists  in  fairly  long  lived  structures, 
most  of  which  is  in  the  form  of  coronal  loops.  The  observed  stability  of  these  struc¬ 
tures,  as  well  as  the  low  0  character  of  the  corona,  suggests  that  they  are  force-free 
configurations,  :.e.,  configurations  in  which  the  magnetic  field.  B.  satisfies  the  equa¬ 
tion  V  X  B  =  oB,  with  a  =  a(x.  t).  Hence  the  disturbances  in  the  coronal  magnetic 
field  induced  by  the  motions  of  the  footpoints  can  be  thought  of  as  perturbations  on 
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a  force-free  state.  The  relative  stability  of  the  coronal  magnetic  field  suggests  that 
the  perturbed  system  then  relaxes  to  a  new,  force-free  state.  Noting  this,  Heyvaerts 
and  Priest  (1984)  adapted  to  the  solar  corona  a  theory  developed  by  Taylor  (1974) 
for  relaxation  to  force-free,  minimum  energy  states  in  the  reversed  field  pinch,  a 
magnetic  confinement  fusion  device.  The  details  of  Taylor’s  theory  are  provided  in 
their  paper  (see  also  the  recent  review  by  Taylor  1986).  Assuming  that  the  coronal 
magnetic  field  is  helical  and  approximately  force-free,  Heyvaerts  and  Priest  used 
Taylor’s  method  to  determine  to  what  state  the  perturbed  coronal  magnetic  field 
would  releoc.  Using  the  assumption  that  the  kinetic  energy  and  internal  energy  are 
negligible  with  respect  to  the  magnetic  energy,  they  minimized  the  total  energy  in 
the  magnetofluid  subject  to  the  constraint  that  the  magnetic  helicity  remain  con¬ 
stant,  and  obtained  constant  a,  force-free  minimum  energy  states.  The  magnetic 
energy  dissipated  zis  heat  could  then  be  determined  by  comparing  the  initial  mag¬ 
netic  field  with  the  minimum  energy  state  magnetic  field.  The  heating  is  produced 
by  the  dissipation  of  magnetic  energy  at  electric  current  sheets  which  are  formed 
locally  in  the  corona.  These  ide2is  were  extended  and  developed  by  Browning  et  al. 
(1986)  and  Browning  amd  Priest  (1986). 

The  theories  of  Taylor  (1974)  and  Heyvaerts  and  Priest  (1984)  are  based  on 
the  stability  of  the  final  magnetic  field  configuration,  with  no  consideration  of  the 
evolution  of  the  magnetofluid  toward  a  force-free  state.  This  dynamic  relaxation 
process  is  now  thought  to  occur  eis  a  consequence  of  the  self-organization  of  the 
helically  turbulent  magnetofluid.  Because  our  results  are  interpreted  in  terms  of  this 
relaocation  theory,  we  review  certain  relevant  details  in  Sect.  2  (see  also  the  recent 
review  by  Hasegawa  1985  and  references  therein).  We  note  here  certain  artificial 
constraints  which  the  Taylor- Heyvaerts  theory  places  on  the  magnetofluid.  This 


Taylor-Heyvaerts  anzdysis  assumes  that  the  kinetic  energy  will  remain  very  small 
with  respect  to  the  magnetic  energy  at  cill  times  of  the  magnetofluid  evolution, 
and  that  a  will  tend  toward  constant  values.  However,  Low  (1982)  has  suggested 
that  non-constant  a,  force-free  magnetic  fields  are  more  common.  This  follows  as  a 
consequence  of  the  fact  that  the  constant  a  magnetic  field  obeys  the  linear  vector 
Helmholtz  equation  V^B  -b  =  0,  which  has  a  smaller  set  of  solutions  that  is 
contzdned  in  the  set  of  solutions  of  the  more  general  force-free  field  equation,  V  x  B  = 
aB,  where  a  can  vary.  Furthermore,  J.  Dahlburg  et  al.  (19S6i,  1987)  report 
the  observation  of  significant  kinetic  energy  formation,  followed  by  the  attainment 
of  a  force-free,  lower  energy  state  characterized  by  nonconstant  a,  in  numerical 
simulations  of  the  reversed-field  pinch  fusion  device.  Hence,  it  is  important  to 
follow  the  nonlinear  evolution  of  a  fully  three-dimensional,  unbounded,  corona-like 
magnetofluid  to  determine:  (1)  if  such  magnetofluids  tend  to  evolve  through  a 
succession  of  force-free  states  when  perturbed  away  from  a  force-free  state;  (2)  if 
the  approach  to  such  states  results  in  enhanced  heat  production;  (3)  if  the  velocity 
field  plays  a  significant  role  in  the  dissipative  process;  and  (4)  if  any  lower  energy 
states  which  occur  are  also  constant  a  states.  We  also  want  to  observe  the  structural 
consequences  of  the  nonlinear  evolution  of  the  magnetofluid,  and  to  determine  what 
structures  contribute  to  the  heating. 

We  here  report  results  from  fully  three-dimensional,  single-fluid  MHD  simula¬ 
tions  in  periodic  geometry  which  help  to  answer  these  questions.  We  consider  an 
initial  value  problem  with  a  helical  force-free  magnetic  field  perturbed  by  random 
noise  in  the  velocity  and  magnetic  fields.  Study  of  the  initial  value  problem  is  sim¬ 
ilar  in  spirit  to  the  "mixing  time"  approach  of  Heyvaerts  and  Priest  (  10S4).  The 
system  is  assumed  to  be  periodic,  so  no  effects  due  to  confinement  are  possiljle. 


We  follow  the  evolution  of  this  system  and  determine  the  extent  and  character  of 
the  accelerated  heating  which  occurs  due  to  turbulent  transfer  of  excitations.  In 
Sect.  3  we  give  the  initial  conditions  and  boundary  conditions.  The  results  of  the 
numerical  simulations  ^lre  described  in  Sect.  4.  In  Sect.  5  we  discuss  the  results. 
Our  numerical  method  is  described  in  the  appendix. 


2.  Review  of  turbulent  self-organization 


The  nonlinear  partial  differential  equations  that  govern  the  behavior  of  a  three- 
dimensionJil,  viscoresistive,  incompressible  magnetofluid,  written  in  a  dimensionless, 
rotation  form,  eire: 


^  =  vxu;-Vn  +  JxB  +  i^V2v  =  H, 
at  M 


(la) 


V-v  =  0, 


(16) 


Q\  1 

^=vxB  +  V$  +  ^V2a  =  G,  (Ic) 

at  b 

V-A  =  0.  (Id) 

where  v(x,  t)  =  flow  velocity,  u)(x,t)  =  V  x  v  =  vorticity,  A(x,  t)  =  magnetic 
vector  potentiad,  B(x,  t)  =  V  x  A  =  magnetic  field,  J(x,t)  =  V  x  B  =  electric 
current  density,  II(x,  t)  =  pressure  head,  $(x,  t)  =  a  magnetic  scalar  potential. 
S  =  Lundquist  number,  and  M  =  viscous  Lundquist  number.  The  Coulomb  gauge 
is  aissumed.  In  this  representation,  B  is  measured  in  terms  of  a  characteristic 
field  strength,  e.g.,  Bq  =  <|  B  j"  ’>^  ,  where  the  brackets  <>  indicate  an 

integral  over  the  system  and  V  denotes  the  volume  of  integration.  The  velocities 
Eue  measured  in  units  of  the  Alfven  speed,  Ca  =  — where  tlie  ma.s.s  density, 
p,  is  assumed  to  be  constant  and  uniform.  For  a  characteristic  distance.  Lo.  the 
time  is  measured  in  units  of  the  .\lfven  transit  time,  The  Lundqtiist  number 


S  =  ,  where  77  is  the  magnetic  difFusivity.  The  viscous  Lundquist  number 

M  =  ,  where  v  is  the  kinematic  viscosity.  The  Lundquist  numbers  express  the 

ratio  between  the  convective  terms  and  the  dissipative  terms  (for  a  fuller  discussion 
of  the  importeince  of  the  Lundquist  numbers,  see  R..  Dahlburg  tt  al.  (1983)). 

For  the  idezd  magnetofluid,  S  —*  oo  and  M  —*  00.  In  this  limit  there  is  no 
dissipation.  There  are  known  to  be  three  conserved,  global  quantities  for  the  three- 
dimensional,  ideal  magnetofluid  {e.g.,  Frisch  ei  al.  1975).  These  are  the  total 
energy: 

2ir  2ir  2jr 

=  J  +  \B\')  dxdydz,  (2a) 

000 


the  magnetic  helicity: 


2ir  2-n 

0  0  0 


B  dx  dy  dz. 


{2b) 


and  the  cross  helicity  (sometimes  called  the  cross  energy): 


Kc  =  ///v 


B  dx  dy  dz. 


-"c 


where  periodic  boundary  conditions  are  cissumed.  These  quantities  also  remain 
invariant  if  perfectly  conducting  and/or  free-slip  boundary  conditions  are  assumed 
in  any  or  all  of  the  three  directions. 

These  ideal  constants  of  motion  are  of  great  importance  m  the  statistical  theory 


of  MHD  turbulence  because  their  Fourier  transforms  define  hypersurfaces  in  the 
complex  Fourier  space  of  coefficients  of  the  magnetic  and  %’elocity  fields.  The  volume 
averaged.  Fourier  representations  of  the  ideal  invariants  are; 


(3a) 


“  k  ^ 

ifM  =  X^[A-(k,0-B(k,<)],  (36) 

k 

and, 

^fc  =  ^[v(k,/)-B(k,0],  (3c) 

k  ^ 

where  the  superscript  denotes  the  complex  conjugate  and  k  is  the  Fourier 
wavenumber. 

The  behavior  of  the  turbulent  ideal  magnetofluid  will  depend  strongly  on  the 
relative  magnitudes  of  these  invariants.  For  e.xample,  consider  the  case  in  which  the 
magnetic  energy  is  much  greater  than  the  kinetic  energy,  the  situation  in  the  solar 
corona.  This  implies  that  the  cross  helicity  will  be  negligible,  since  the  magnitude 
of  the  velocity  field  will  be  small.  Furthermore,  the  total  energy  will  consist  almost 
exclusively  of  magnetic  energy,  since  in  this  ceise  |  B  |^>>|  v  |^. 

In  the  presence  of  finite  dissipation,  the  total  energy,  Er.  and  the  magnetic 
helicity,  Hm,  will  decay  according  to  the  following  equations  (e.g..  Matthaeus  and 
Montgomery  1980): 

^ = -1 7 77 1 J  r  -  A  7// 1 .  r 

0  0  0  0  0  0 

2Uld, 


7 


(46) 


^=-i/  1 1  IBd^dydz. 

0  0  0 

The  first  integral  on  the  right  hand  side  of  equation  4a  is  called  the  mean  square 
electric  current,  and  it  determines  the  rate  at  which  the  magnetic  energy  is  dis¬ 
sipated.  The  second  integral  on  the  right  hand  side  of  equation  3a  is  called  the 
enstrophy,  amd  it  determines  the  rate  at  which  the  kinetic  energy  is  dissipated.  For 
the  I  B  |^>>|  V  1^  Ccise,  the  enstrophy  is  much  smailler  than  the  mean  square  elec¬ 
tric  current.  Hence,  for  a  unit  magnetic  Praindtl  number  the  viscous  dissipation  of 
energy  should  be  negligible.  The  Fourier  space  representation  of  these  decay  laws 
is: 


k  k 

and, 

^  =  (56) 

k 

Now  consider,  for  example,  a  force-free  case  in  which  the  magnetic  vector  poten¬ 
tial,  the  magnetic  field,  and  the  electric  current  density  are  all  directly  proportional 
to  each  other  initially  (we  will  consider  such  a  case  in  this  paper).  In  this  case  the 
ratio  Et/H^{  will  remain  constant  as  the  magnetic  field  decays  by  means  of  Ohmic 
diffusion.  When  perturbed,  however,  the  system  will  execute  a  more  complex  evo¬ 
lution  due  to  the  activation  of  the  nonlinear  terms  in  the  governing  (Mpiations.  As 
the  system  evolves  nonlinearly  the  magnetic  field  excitations  will  tc'nd  to  sjna'ad 


out  in  wavenumber  space  due  to  interactions  among  resonant  triads  of  wavevectors. 
The  electric  current  density,  J(k)  will  populate  the  shorter  wavelengths  at  a  rate 
fjister  them  the  magnetic  field,  B(k),  since  J{k)  =tk  x  B(k).  Fort  >  l,J*(k)- J(k) 
will  tend  to  be  larger  in  magnitude  than  J*(k)  ■  B(k).  Furthermore,  J  •  J  is  posi¬ 
tive  definite,  while  J  •  B  can  be  negative.  All  this  implies  that  a  selective  decay  of 
the  magnetic  energy  will  take  place  with  respect  to  the  magnetic  helicity  once  the 
turbulent  transport  of  excitations  becomes  significant  (Montgomery  et  al.  197S). 
Consequently  the  ratio  Et/Hm  should  decreeise  as  a  function  of  time  when  the 
system  becomes  turbulent. 

Finzdly,  the  nonlinear  mode  coupling  is  such  that  the  magnetic  helicity,  which 
suffers  less  dissipation  than  the  total  energy,  is  preferentially  transferred  to  the 
longest  wavelengths  of  the  system  (Frisch  et  al.  1975,  Pouquet  et  al.  1976.  Mont¬ 
gomery  et  al.  1978).  This  transfer  of  magnetic  helicity  to  smaller  wavenumbers  finds 
morphological  expression  in  the  tendency  of  the  magnetic  field  to  assume  long- 
wavelength  structures.  The  transfer  of  total  energy  to  larger  wavenumbers  finds 
morphological  expression  in  the  formation  of  electric  current  sheets  and  concen¬ 
trated  vorticity  structures.  These  current  sheets  and  concentrated  vorticity  struc¬ 
tures  are  the  primary  sites  at  which  turbulent  heating  occurs. 

The  finzd  state  the  relaxing  magnetofiuid  is  thought  to  assume  is  a  minimum- 
energy,  force-free  state,  commonly  referred  to  as  a  "Taylor  state."  These  final  states 
have  been  discussed  by  Taylor  ( 1974)  and  Heyvaerts  and  Priest  ( 19S4),  who  consid¬ 
ered  stability  properties  alone.  By  minimizing  the  magnetic  energy  su!:>ject  to  the 
constraint  of  constant  magnetic  helicity.  Taylor  showed  that  the  minimum  energy 
state  to  be  obtained  would  be  a  constant  a  state,  where  V  x  B  =  oB.  .\lthough  the 
evolution  toward  a  force-free  sttrte  has  been  found  to  occur  in  computer  simulations 


by  J.  DaMburg  et  ai  (1987)  in  reversed-field  pinch-like  geometry,  questions  remain 
with  respect  to  the  low  velocity  and  constzmt  a  assumptions,  zis  well  as  the  possi¬ 
bility  of  evolution  toward  a  force-free  state  in  situations  in  which  rigid,  perfectly 
conducting,  impenetrable  walls  are  not  present,  e.g.,  in  the  solau:  corona. 


3.  Formulation  of  the  problem 


For  the  initizil  background,  force-free  magnetic  field  we  use  the  following  pre¬ 
scription  of  the  magnetic  vector  potential: 

-4i(x,  f  =  0)  =  Cl  sinicosy  ,  (6a) 

Ay(x,t  =  0)  = —Cl  cos  X  siny  ,  (66) 

Aj(x,  t  =  0)  =  C2  sinx  siny  ,  (6c) 

where  Ci  =  0.5  and  C2  =  C°'^.  No  current  sheets  are  present  initially  in  this 
magnetic  field  configuration.  Note  that  initially  V  x  B  =  oB,  with  a  =  2“  ^ 
throughout  the  system  (  in  fact,  V'B  -I-  a^B  =  0,  and  this  is  a  linear,  constant 
a,  force-free  magnetic  field).  The  initial,  unperturbed,  volume  averaged  magnetic 
energy  equals  1.  The  initial,  unperturbed,  volume  averaged  magnetic  helicity  equals 
i  _  9-0.5 

a  ^ 

This  configuration  (6a  —  6c)  will  simply  decay  Ohmically  if  unperturbed,  since 
fluid  motion  will  be  excited  because  it  is  force-free.  We  perturb  this  configuration 
with  broad-band  noise  in  both  the  velocity  and  magnetic  field  to  ensure  a  nontrivial 
evolution.  This  allows  the  system  to  traverse  as  wide  a  range  of  the  availaide  pliase 
space  as  possible.  The  fastest  growing  disturbances,  whatever  they  might  lie.  emerge 
naturailly  from  the  noise  as  time  advances.  These  initi.al  conditions  are  similar  to 
those  used  by  Horiuchi  and  Sato  ( lOSo.  lOSG). 


Periodic  boundary  conditions  are  enforced  in  all  three  spatial  directions.  It 
seems  wisest  to  use  these  boundziry  conditions  so  that  the  turbulent  magnetofluid 
can  be  simulated  with  the  greatest  accuracy,  i.e.,  without  the  need  to  use  special 
measures  to  enforce  more  specific  boundary  conditions  {cf.  Meneguzzi  et  al.  1981, 
Pouquet  et  al.  1986).  The  restriction  to  periodic  boundary  conditions,  coupled 
with  our  use  of  a  Fourier  pseudospectral  numerical  algorithm,  allow  us  to  obtain 
high  resolution  of  the  turbulent  evolution  of  the  magnetofluid.  The  fully  three- 
dimensional,  Fourier  pseudospectral,  numerical  simulation  code  used  to  advance 
the  governing  equations,  la  —  Id,  is  described  in  the  appendix. 


4.  Results  of  numerical  simulations 


In  this  section  we  describe  the  results  from  a  typical  simulation,  which  is  rep¬ 
resentative  of  about  ten  simulations  that  we  have  performed.  The  spatial  grid  uses 
32*  Fourier  modes,  with  a  time  step  of  1/500.  At  this  size  the  code  requires  about 
1.5  seconds  per  time  step  on  the  CRAY-2,  and  approximately  1.2  million  words 
of  memory.  The  Lundquist  numbers  are  set  to  400  for  this  run.  On  the  baisis  of 
the  Kolmogorov  dissipation  wavenumber,  this  is  close  to  the  limit  at  which  we  can 
accurately  compute  the  turbulent  evolution  of  the  magnetofluid.  In  fact,  accurate 
numerical  simulation  of  two-dimensional  turbulent  magnetofluids  hzis  been  limited 
to  Lundquist  numbers  of  about  1000  on  the  present  generation  of  computers  (e.g., 
Forbes  and  Priest  1983)  and  to  Lundquist  numbers  of  less  thEin  500  for  the  three- 
dimensional  case  (J.  Dalilburg  et  al.^  1987).  These  Lundquist  numbers  are  well 
above  the  known  linear  stability  boundaries  for  magnetic  reconnection  to  occur  (R. 
Dahlburg  et  al.  1983,  Bondeson  and  Sobel  1984). 

We  first  describe  the  evolution  of  the  magnetofluid  with  the  assistance  of 
the  global  diagnostics,  and  then  discuss  the  structural  evolution.  Figure  1  shows 
Eb,Ev,Hm,  and  Et/Hm  as  functions  of  time.  All  global  quantities  shown  are 
volume  averaged.  The  magnetic  energy  in  the  configuration  described  by  equations 
6a  -  6c  is  equzd  to  1.  The  initial  kinetic  energy  and  the  perturbed  magnetic  energy 
each  are  equal  to  0.25x10“*,  with  the  band  of  excitation  given  by  9  <|  k  |‘<  1000. 
Hence  the  energy  in  the  perturbation  is  0.05  %  of  the  background  magnetic  configu¬ 
ration  given  by  equations  6a  —  6c.  The  magnetic  energy,  Eb^  first  decays  Ohmically 
for  about  twenty  Alfven  transit  times.  At  this  time,  turbulent  activity  begins  and 
accelerated  dissipation  of  Eb  occurs.  The  lost  magnetic  energy  is  both  dissipated 
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and  transformed  into  kinetic  energy.  Following  this,  some  dynamo  activity  is  ob¬ 
served  as  the  magnetic  energy  resumes  decaying  at  a  slower  rate.  The  kinetic  energy, 
Evi  remains  relatively  small  until  the  accelerated  decay  of  the  magnetic  energy  oc¬ 
curs.  At  this  time  a  rapid  buildup  of  kinetic  energy  occurs,  indicative  of  organized 
fluid  motions.  After  severed  bursts  of  activity,  the  kinetic  energy  dies  off.  The 
magnetic  helicity,  Hm,  decays  relatively  slowly  throughout  the  run,  with  relatively 
little  variation  in  response  to  the  magnetofluid  activity.  Of  perhaps  more  impor¬ 
tance  is  the  ratio  of  the  toted  energy  to  the  magnetic  helicity,  In  figure 

1  we  normedize  this  ratio  to  its  initial  vedue  of  1.414.  This  ratio  remains  appro.xi- 
mately  constant  until  the  onset  of  turbulent  activity,  indicating  that  the  decay  of 
the  toted  energy  emd  the  magnetic  helicity  at  eeurly  times  is  being  determined  by  the 
lineeir,  dissipative  terms.  The  decline  in  the  value  of  this  ratio  indicates  the  onset 
of  turbident  activity,  in  which  the  total  energy  is  more  effectively  dissipated  than 
the  magnetic  helicity.  This  occurs  because  energy  is  being  transferred  to  shorter 
wavelengths  where  dissipation  is  more  effective,  and  because  am  additional  energy 
dissipation  process,  viscous  dissipation,  is  now  available  to  the  magnetofluid.  .A.t 
later  times  this  ratio,  Et/Hji^i,  appears  to  asymptote  toward  a  lower  value  than  it 
had  initially. 

More  exact  information  about  the  energy  dissipation  is  shown  in  figure  2,  which 
shows  the  mean  square  electric  current,  the  enstrophy,  and  the  fractional  heating 
rate  as  functions  of  time.  All  three  of  these  quantities  are  seen  to  rise  rapidly  at  the 
onset  of  the  turbulent  heating  process.  The  electric  current  density  and  the  vorticity 
empheisize  the  high  k  parts  of  the  magnetic  field  and  velocity  field,  respectively. 
Hence  the  mean  square  electric  current  and  the  enstrophy  grow  as  the  magnetofluid 
becomes  turbulent.  We  will  show  later  that  the  rapid  rise  in  the  mean  square  electric 
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current  is  coincident  with  the  formation  of  sheets  of  electric  current.  In  the  same 
way,  the  rise  in  the  enstrophy  is  coincident  with  the  generation  of  concentrated 
vorticity  structures  in  the'  vicinity  of  these  electric  current  sheets.  The  fractional 
heating  rate,  i.e.,  the  logarithmic  derivative  of  the  toted  energy,  equals  the  sum  of 
the  enstrophy  and  mean  square  electric  current  divided  by  the  total  energy.  This 
rate  remains  approximately  constant  until  the  onset  of  turbulence,  indicating  that 
the  dissipation  up  to  this  time  is  not  being  mediated  by  nonlinear  processes.  When 
the  turbulence  sets  in,  a  rapid  rise  in  the  fractionzd  heating  rate  occurs,  followed  by 
several  smaller  bursts  of  activity. 

It  is  of  interest  that  the  enstrophy  attains  such  a  high  value  relative  to  the  mean 
square  electric  current.  Figure  3  shows  the  ratio  of  kinetic  energy  to  magnetic 
energy  as  a  function  of  time,  and  the  ratio  of  the  enstrophy  to  the  mean  square 
electric  current.  The  initial  value  of  the  ratio  Ev/Eb  is  2.5x10—*,  and  it  attains 
a  maximum  value  of  about  0.1.  The  ratio  of  the  viscous  dissipation  to  the  Ohmic 
dissipation  has  an  initial  value  of  about  1.0x10“*  and  attains  its  maximum  value 
of  almost  0.7.  Both  of  these  ratios  attain  their  peak  values  at  about  the  time  of 
greatest  heating  {t  w  30)  ais  determined  by  the  fractional  heating  rate  of  figure 
2.  Hence,  although  the  kinetic  energy  is  small,  the  turbulence  of  the  magnetofiuid 
produces  a  relatively  lauge  amount  of  viscous  dissipation.  This  also  implies  that 
a  significant  zunount  of  heating  occurs  by  a  secondary  process,  viz.,  the  magnetic 
energy  is  transformed  first  into  kinetic  energy  and  then  into  heat  energy  by  viscous 
dissipation. 

Figure  4  exhibits  the  harmonic  history  for  several  of  the  magnetic  energy  modes. 
Initially,  most  of  the  magnetic  energy  is  in  the  (1,1.0)  mode  {i.e...  (kj.  =  l.k,^  = 
l,ki  =  0)),  for  which  |  k  |’=  2.  The  initial  ])erturbation  of  the  magnetic  field  is 
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entirely  in  modes  with  |  k  |  >2,  and  no  energy  is  in  any  mode  with  |  k  |  <2.  The 
(1,1,0)  mode  decays  Ohmically  until  the  onset  of  MHD  turbulence.  After  about 
twenty  Alfven  transit  times  a  period  of  rapid  decay  of  the  (1,1,0)  magnetic  mode 
occurs.  At  this  time  the  energy  is  being  transferred  out  of  the  (1,1,0)  mode  to  modes 
with  different  values  of  |  k  in  both  the  magnetic  field  cind  the  velocity  field  by 
turbulent  processes.  At  the  shorter  wavelengths  the  dissipation  is  more  effective. 
We  tilso  see  the  excitation  of  certain  longer  wavelength  modes  of  the  system,  in 
particular  the  (1,0,0)  mode  and  the  (0,1,0)  mode.  For  both  of  these  modes  [  k  |^=  1. 
Within  the  period  of  a  few  Alfven  transit  times  these  modes  grow  so  rapidly  that 
their  energies  each  exceed  in  magnitude  the  (1,1,0)  mode.  Excitation  of  the  (0,0,1) 
mode  is  comparatively  negligible.  The  morphological  consequence  of  this  turbulent 
transfer  of  energy  is  shown  in  figure  5,  which  shows  contours  of  constant  z  magnetic 
field  at  times  early  zmd  late  in  the  numerical  simulation.  This  e.xcitation  of  the 
lower  I  k  modes  is  similar  to  that  derived  by  Browning  et  al.  (19S6). 

The  processes  which  contribute  to  the  accelerated  heating  eure  illuminated  by 
examination  of  two-dimensional  slices  taken  through  the  system  perpendicular  to 
the  z  direction  at  (x,y,z  =  0).  For  several  important  times  we  show  plots  of  Di 
and  By,  7*,  u  amd  u,  and  w..  Figure  6  shows  these  fields  at  ^  =  0.5,  a  time  prior  to 
the  onset  of  the  rapid  heating.  The  state  of  the  fields  at  this  time  is  very  similar  to 
the  initiad  conditions.  The  magnetic  field  By)  and  electric  current  density  ( J; ) 
are  similar  in  appearance  due  to  the  almost  force- free  character  of  the  magnetic 
field.  However,  the  electric  current  density  (7;)  appears  more  perturl>ed  Ijecause 
it  emphasizes  the  high  wavelength  components  of  the  magnetic  field.  These  fields 
cire  close  to  zero  in  the  centre  of  the  plot  and  at  the  corners.  No  electric  current 
sheets  are  present  initially.  The  velocitv  field  (a,  v)  and  vorticitv  (u-'A  are  relativelv 


random  in  appe2irance.  This  also  reflects  their  initial  conditions.  Note  the  absence 
of  long-wavelength  structure  in  the  velocity  fleld  (u,  y)  and  vorticity  field  {uj; ),  and 
the  lack  of  correlation  with  the  magnetic  field  (Br,By). 

Figure  7  shows  these  fields  at  t  =  25,  during  the  phcise  of  rapid  heating.  Two 
changes  in  the  magnetic  field  (Bz,  By)  configuration  axe  apparent,  viz.,  the  shifting 
of  the  field  so  that  regions  of  oppositely  directed  field  are  adjacent,  and  emergence 
of  larger  scales  in  the  magnetic  field.  One  merger  is  occuring  at  the  center  of  the 
plots,  the  other  is  occuring  at  the  corners.  Since  these  mergers  are  similar  we  will 
restrict  discussion  to  the  one  occuring  in  the  centre  of  the  plot.  Note  that  since  the 
magnetic  field  is  not  confined  by  walls  it  is  free  to  reconfigure  itself  in  this  fashion. 
The  merging  of  the  magnetic  field  is  mediated  by  the  large  electric  current  sheet 
which  forms  between  the  merging  regions.  This  sheet  ceirries  electric  current  of  sign 
opposite  to  the  electric  current  in  the  merging  regions.  The  maximum  amplitude  of 
the  electric  current  density  (J,)  in  the  sheet  is  four  to  five  times  that  of  the  initial 
maxima.  It  is  the  formation  of  these  sheets  that  leads  to  the  increase  in  the  mean 
squeire  electric  current  at  this  time  (figure  2).  Large  scale,  organized  fluid  motions 
also  are  now  seen  in  (u,y).  These  fluid  motions  serve  to  drag  magnetic  flux  into 
the  electric  current  sheet,  increasing  the  reconnection  rate  and  hence  the  heating 
rate.  Note  eilso  the  flows  out  of  the  reconnection  zone,  with  speeds  a  significant 
fraction  of  the  Alfven  speed.  A  quadrupolar  vorticity  distribution  in  (a,’;)  is  seen  to 
form  in  the  vicinity  of  the  electric  current  sheet.  The  vorticity  (u,’;)  maxima  at  this 
time  are  over  thirty  times  greater  than  their  initial  value.  The  viscous  dissipation 
of  energy  is  very  strong  in  these  regions  of  large  vorticity.  and  this  is  evidenced  by 
the  increase  in  enstrophy  seen  after  about  t  =20  (figure  2). 


Figure  8  shows  these  fields  at  the  final  time  in  the  simulation,  /  =  53.  At  this 
time  most  of  the  turbulent  activity  has  ceased.  The  magnetic  field  {Bz,By)  is  now 
dominated  by  the  largest  scales  in  the  system.  This  also  can  be  inferred  from  a 
comparison  of  figure  1  zmd  figure  4.  At  this  time  the  electric  current  density  ( Jt ) 
has  decreased  in  value,  the  maccima  being  about  one  third  of  their  original  value. 
Although  the  electric  current  density  is  noisy,  its  form  appears  to  be  similar  to 
that  of  the  magnetic  field  -  a  characteristic  to  be  expected  in  a  force-free  state. 
The  electric  current  sheets  which  produced  the  accelerated  dissipation  have  disap¬ 
peared.  The  velocity  field  (u,u)  still  possesses  some  long  wavelength  structure,  but 
its  magnitude  has  declined.  The  vorticity  field  (u;^ )  now  has  an  appearance  similar 
to  the  electric  current  density  with  no  concentrated  vorticity  structures  evident, 
and  the  maxima  of  both  Jz  and  Uz  are  of  similar  magnitude. 

One  feature  of  interest  is  the  behavior  of  a  in  the  force-free  regions.  We  show 
slice  plots  of  a  and  the  force-free  regions  in  the  (i,y,  z  =  0)  plane  which  indicate 
that  o  varies  significantly  in  the  force-free  regions.  Figure  9  shows  a  and  the  force 
free  regions  at  <  =  0.5  amd  t  =  53.  To  delineate  the  force  free  regions,  we  plot 
(j-E)/(i  j  n  B  p  )K  This  quantity  will  have  a  magnitude  of  1  in  the  force-free 
regions.  To  show  the  spatial  distribution  of  a,  we  plot  (J  •  B)/(B  •  B).  Recall  that 
a  =  2^  for  the  magnetic  field  given  by  equations  6a  —  6c.  .A.t  t  =  0.5  the  magnetic 
field  is  almost  entirely  force  free  except  for  small  regions  near  the  centre  and  at  the 
corners  of  the  plot.  These  are  the  regions  in  which  B  is  close  to  zero.  At  the  same 
time,  some  variation  in  a  is  seen  throughout  the  force-free  regions.  The  variation  in 
a  is  due  to  the  initial  perturbation,  and,  as  might  be  expected,  the  largest  variation 
is  seen  in  the  regions  which  are  not  force-free.  .A.t  t  =  53,  the  magnetic  field  has 


evolved  to  a  new  state  which  has  large  regions  that  axe  force-free.  Substantially 
more  variation  in  a  is  evident,  especially  near  the  corners  of  the  plot. 

We  have  performed  these  simulations  at  S  =  M  =  50, 100,  cind  200  with 
essenti2dly  the  same  results.  In  the  absence  of  nonlinear  effects,  rausing  S  slows  down 
the  decay  of  the  background  magnetic  field.  In  the  presence  of  nonlinear  effects, 
raising  the  Lundquist  numbers,  M  and  S  increases  the  strength  of  the  nonlinear 
terms  relative  to  the  dissipative  terms.  This  edlows  the  magnetofluid  to  t.-ansfer 
excitations  to  short  wavelengths  more  efficiently.  Hence,  as  the  Lundquist  numbers 
are  raised,  the  peak  magnitudes  of  the  electric  current  sheet  and  the  concentrated 
vorticity  structures  tend  to  increase.  This  leads  to  lairger  peak  values  of  the  mean 
square  electric  current  and  enstrophy,  with  a  corresponding  larger  loss  of  magnetic 
energy  during  the  phase  of  rapid  dissipation.  The  heating  contribution  due  to 
enstrophy  generation  is  significant  in  all  of  these  runs.  Variation  in  a  also  was 
seen  throughout  all  of  the  runs  performed.  We  eiIso  have  used  other  sets  of  random 
perturbations,  and  ageiin  obtained  essentially  the  same  results.  A  full  parameter 
search  is  impractical  because  of  the  expense  of  these  calculations. 


5.  Discussion 


In  this  paper  we  report  the  results  of  numerical  simulations  of  the  decay  of 
ttirbulence  in  a  maguetofluid  with  a  high  degree  of  magnetic  helicity  present  ini¬ 
tially.  Similar  magnetofluid  configurations  are  believed  to  be  present  in  the  solar 
corona,  aind  perhaps  to  underlie  the  process  by  which  the  solar  corona  is  heated 
(c/.  Heyvaerts  and  Priest  19S4).  Our  motivation  for  performing  these  numerical 
simulations  was  to  understand  better  the  dynamical  evolution  and  energetics  of  a 
heliceJly  turbulent  magnetofluid,  and  to  identify  the  structures  which  develop  dur¬ 
ing  the  course  of  this  evolution.  For  the  foreseeable  future,  numerical  simulations 
appear  to  be  the  best  way  of  obtaining  such  insight,  since  direct  measurement  of  the 
magnetic  field  and  velocity  field  in  the  solar  corona  is  not  possible  at  present.  The 
ideal  invariants  (equations  2a  —  2c)  have  been  mezisured  in  the  solar  wind  by  means 
of  in  situ,  single  point  covariances  {e.g.,  Matthaeus  and  Goldstein  19S2).  However, 
the  most  common  methods  for  determining  the  coronal  magnetic  field  and  velocity 
field  are  indirect.  The  coronal  magnetic  field  is  commonly  extrapolated  from  the 
longitudinal  component  of  the  photospheric  field,  under  the  assumption  that  it  is 
either  a  potentied  magnetic  field  or  a  force-free  magnetic  field.  The  coronal  velocity 
field  generally  is  extrapolated  from  Doppler-shift  measurements  of  the  line-of-sight 
velocity  field  component  in  the  transition  region.  Simultaneous  specification  of 
both  of  these  fields  at  the  same  point  in  the  corona  is  at  present  out  of  the  ques¬ 
tion.  Analysis  of  helically  turbulent  magnetofluids  is  also  difficult,  primarily  due  to 
the  complexity  of  the  governing  nonlinear,  partial  differential  equations.  Even  very 
simple,  deterministic,  initial  conditions  can  lead  to  analytically  intractable  behav¬ 
ior  after  a  few  Alfven  transit  times  {e.g..  Pouquet  et  al.  19S6).  promising  start 


in  the  analytic  direction  hats  been  made  by  Krishan  (1985),  who  hais  formulated  a 
statistical  mech^lnics  for  active  region  coroneil  magnetofluids.  In  light  of  the  state  of 
observations  and  analysis,  we  have  chosen  not  to  concentrate  our  efforts  on  the  in¬ 
vestigation  of  em  elaborate  model.  Instead,  we  have  investigated  a  relatively  simple 
magnetohydrodynamic  configuration  which  possesses  the  most  significant  feature. 
viz.,  twisted  magnetic  fields,  and  have  attempted  to  understand  the  important  fea¬ 
tures  of  its  evolution.  This  approach  concedes  two  present-day  limitations  of  direct 
numerical  simulations  of  turbulent  fluids;  the  restriction  to  simple  geometries  and 
the  restriction  to  low  Lundquist  (or  Reynolds)  numbers. 

Our  numericcil  simulations  have  shown  that  accelerated  dissipation  of  magnetic 
energy  occurs  in  a  magnetofluid  with  a  high  degree  of  magnetic  helicity  present 
initially  and  with  open  boundary  conditions.  Large  regions  remain  force-free  during 
the  entire  evolution  of  the  magnetofluid.  The  initial  magnetic  configuration  has  no 
electric  current  sheets  present,  but  it  evolves  toward  a  state  in  which  electric  current 
sheets  do  exist.  In  this  respect,  the  boundary  conditions  seem  to  be  important 
because  they  rdlow  the  “tubes”  to  slide  around  each  other.  Formation  of  the  electric 
current  sheet  produces  non-negligible  Lorentz  forces.  High  speed  flows  develop 
which  sweep  additional  magnetic  flux  into  the  region  of  the  electric  current  sheet, 
leading  to  an  accelerated  dissipation  of  magnetic  energy.  It  is  in  the  vicinity  of 
the  non-force-free  field  regions  that  the  bulk  of  the  accelerated  heating  occurs.  The 
final  form  of  the  magnetic  field  appears  to  be  dominated  by  the  largest  length  scales 
rdlowed  in  the  system  (compare  figure  1  and  figure  4).  The  structural  evolution  of 
the  magnetic  field  and  the  electric  current  density,  comljined  with  the  selective  decay 
of  the  total  energy  with  respect  to  the  magnetic  helicity  (figtire  1),  suggest  that  the 
magnetofluid  is  self-organizing  in  the  manner  ilescriljed  in  ‘^ect.  2. 


We  also  note  that  the  magnetofluid  velocity  plays  a  significant  role  in  the 
evolution  of  the  system.  When  the  Lorentz  force  becomes  non- negligible,  high  speed 
flows  develop.  The  highest  level  of  kinetic  energy  is  achieved  during  the  time  of  the 
formation  of  the  electric  current  sheets,  at  which  time  the  kinetic  energy  is  about  10 
%  of  the  magnetic  energy.  The  high  level  of  kinetic  energy  does  not  seem  to  inhibit 
the  self-organization  of  the  magnetofluid,  e.g.,  the  totcil  energy  decays  selectively 
with  respect  to  the  magnetic  helicity  in  a  monotonic  fashion  (figure  1).  However,  the 
numerical  simulations  imply  that  the  magnetofluid  cannot  be  modelled  as  evolving 
quasi-staticzflly  (e.g.,  Heyvaerts  and  Priest  19S4),  and  that  fluid  motions  must  be 
taken  into  account. 

Accompanying  the  increase  in  kinetic  energy  is  an  increase  in  enstrophy,  which 
peaks  at  about  70  %  of  the  mean  square  electric  current  during  the  time  of  greatest 
heating.  Hence,  viscous  heating  aJso  contributes  significantly  to  the  accelerated  dis¬ 
sipation  rate  (c/.  Rosner  et  al.  1986).  We  have  performed  our  numerical  simulation 
with  unit  magnetic  Prandtl  number.  However,  in  the  solar  corona  the  magnitude  of 
the  kinematic  viscosity  is  estimated  to  exceed  the  magnetic  diffusivity  {e.g..  Pneu- 
man  amd  Orrall  1986).  This  implies  that  the  heating  due  to  viscous  dissipation 
might  be  greater  than  that  due  to  Ohmic  dissipation.  The  rise  in  enstrophy  is 
directly  related  to  the  formation  of  concentrated  structures  in  the  v'orticity,  which 
in  this  caise  take  the  form  of  a  “quadrupolar  prism”  (similar  in  form  to  structures 
found  in  two  dimensions,  cf.  Matthaeus  1982).  These  structures  are  similar  to  the 
electric  current  sheets  with  respect  to  dissipation,  viz.,  within  them  th>’  vorticify 
cam  become  so  large  that  even  at  very  low  viscosities  it  is  possible  to  dissip  itv  i...ge 
amounts  of  kinetic  energy.  Our  simulations  appear  to  show  that  these  concentrated 


vorticity  structures  are  cis  important  in  the  coronal  heating  process  as  are  electric 
current  sheets. 

Another  feature  of  interest  is  the  temporal  behavior  of  a,  the  ratio  of  the 
electric  current  density  to  the  magnetic  field.  Our  numerical  results  indicate  that 
a  does  not  approach  a  constant  value  after  the  time  of  rapid  heating,  as  suggested 
by  the  theory  of  Heyvaerts  and  Priest  (1984).  Rather,  although  the  magnetofluid 
quickly  returns  to  a  state  that  can  be  characterized  zis  mostly  force-free,  it  is  not 
observed  to  return  to  a  state  for  which  the  electric  current  density  and  the  magnetic 
field  are  topologically  identical.  Thus,  the  accelerated  heating  occurs  as  the  system 
evolves  through  a  series  of  force-free  states,  but  not  through  a  series  of  constant 
a  states.  We  suggest  that  an  explanation  for  this  may  lie  in  the  direct  action  of 
fluid  motions.  When  the  Lorentz  force  in  the  magnetofluid  becomes  locally  non- 
negligible,  fluid  motions  axe  excited.  Generally,  the  evolving  velocity  field  will  not 
be  aligned  with  the  magnetic  field,  and  hence  it  will  perturb  the  magnetic  field 
away  from  2my  constant  a  force-free  state  through  the  perturbed  electric  field  term. 
V  X  B,  in  the  magnetic  equations.  When  the  Lorentz  force  term  becomes  weaker 
as  a  consequence  of  the  fluid  motions  and  dissipative  losses,  the  fluid  acceleration 
also  decreaises  until  the  J  x  B  force  term  hris  entirely  disappeared  and  the  fluid  is 
again  nearly  at  rest.  The  resulting  nearly  quiescent  state  must  be  force-free  but  not 
necessEirily  constant  a,  which  fits  the  final  states  from  our  simulations.  This  implies 
that  a  constant  a  approximation  must  be  applied  with  great  caution. 

We  are  currently  extending  this  work  to  two  additional  cases  of  interest  to  the 
solar  physics  community.  First,  we  wish  to  determine  the  influence  of  a  constant 
magnetic  field  on  the  results  reported  here.  Second,  we  are  determining  the  evolu¬ 
tion  of  this  system  when  it  is  subjected  to  constant  injection  of  magnetic  helicity. 


both  with  and  without  the  presence  of  a  non-self-consistent  magnetic  field  compo¬ 
nent.  This  situation  is  analogous  to  the  twisting  of  the  coronal  footpoints  by  fluid 
motions  in  the  photosphere  {e.g.,  Parker  1982,  Browning  ei  al.  1986).  Finally,  the 
addition  of  compressible  effects  and  radiation  to  the  model  is  of  interest  (Pettini  et 
al.  1985;  R.  Dahlburg  et  al.  19866,  1987,  Karpen  et  al.  1987). 
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Appendix 


Numerical  algorithm 

Since  we  are  considering  a  fully  periodic  magnetofluid,  it  is  appropriate  to 
use  a  Fourier  pseudospectral  method  for  the  spatial  discretization  (Gottleib  and 
Orszag  1977).  In  the  customary  fashion,  we  evaluate  all  spatiail  derivatives,  without 
phase  error,  in  the  relevant  partial  Fourier  space.  The  nonlinear  product  terms  are 
evaluated  most  easily  in  configuration  space.  We  use  the  vectorized  fast  Fourier 
transform  algorithm  of  Temperton  (1983)  to  move  between  configuration  space  and 
Fourier  space.  We  note  that  pseudospectreil  methods  have  been  applied  with  great 
success  to  a  variety  of  hydrodynamic  and  magnetohydrodynaunic  turbulence  and 
transition  problems  (see  the  recent  review  by  Zang  aind  Hussaini  1987). 

The  equations  with  time  derivatives,  la  and  Ic,  are  written  in  the  rotation 
form.  In  the  ideal  case,  the  discretized  form  of  these  equations  is  said  to  semi¬ 
conserve  energy  -  in  the  sense  that  the  only  errors  in  the  energy  conservation  are 
due  to  the  time  discretization.  A  consequence  of  this  energy  conservation  property 
is  improved  stability  of  the  numerical  method.  A  formal  proof  of  this  point  is  given 
by  T.  A.  Zang  in  an  appendix  of  J.  Daihlburg  et  al.  (19S6a)  (see  also  Gottleib  and 
Orszag  1977). 

For  our  calculation  we  find  it  convenient  to  enforce  the  solenoidality  of  the 
velocity  and  magnetic  vector  potential  fields  indirectly  by  means  of  the  appropriate 
Poisson  equations  (Note  that  the  solenoidality  of  the  magnetic  induction  field  is 
enforced  trivially,  i.e.,  V-B  =  V-  VxA  =  0).  Taking  the  divergence  of  equations 
la  and  Ic  gives  the  result: 


«  -  .  w  .  A  w*.  -V  .  .  ■ 
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.*^j 


v2n  =  VH, 


v2$  =  -v  -  G 


(.-11a) 


(Alb) 


Solving  the  Poisson  equations  Ala  cind  .416  for  11  and  $  ensures  the  solenoidality 
of  the  velocity  cind  magnetic  vector  potential  fields  (see  R.  Dahlburg  et  al.  19S6a). 

In  order  to  explicate  more  fully  the  numerical  method,  we  first  rewrite  equations 
la  and  Ic: 


=  G 


(.-12a) 


(.426) 


We  time  advance  v  and  A  by  a  second  order  Runge  Kutta  scheme.  For  example, 
consider  the  x  component  of  the  velocity  field,  u(x,  t).  For  ^  =  Hi{v.  A,t),  the 


time  discretization  gives: 


=  a"  +  ^/f;(v™,A",t"), 
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(.-136) 


where  the  superscript  n  indexes  the  time  level.  In  the  absence  of  dissipation,  the 
temporal  discretization  is  known  to  be  weakly  unstable,  but  this  should  not  affect 
integration  to  the  time  levels  that  we  consider  here.  A  full  discussion  of  the  stability 
and  accuracy  properties  of  this  time-stepping  scheme  is  given  by  J.  Dahlburg  et  al. 
(1985). 

The  Poisson  equations,  Ala  and  A16,  are  eaisily  inverted  in  the  full  Fourier 
space,  where  they  become: 

fc^n(k,  t)  =  — ik  ■  H(k,  t) ,  (.44a) 

k^^k,t)  =  ik-G{kJ).  (.446) 

In  stimmary,  we  first  solve  for  II  and  $  at  the  previous  time  level  n.  We  then 
time-advance  v  and  A  to  the  time  level  n  +  ^.  Then  we  solve  for  11  and  >5  at  the 
time  level  n  -f  |.  Using  these  values  for  the  time  level  n  -I-  1  we  perform  the  full 
time-step  on  v  and  A  to  the  time  level  a  -I-  1. 
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Figure  1.  Plot  of  global  quantities  us  time.  Dashed  line;  magnetic  energ}':  Dash-<lot 
line:  kinetic  energy  times  ten;  Mixed  dash  line:  magnetic  helicity.  Solid  line: 
ratio  of  toted  energy  to  magnetic  helicity. 


Figure  3.  Plot  of  global  quantities  vs  time.  Solid  line;  Ratio  of  kinetic  energy 
to  magnetic  energy;  Dashed  line;  Ratio  of  enstrophy  to  mean  square  electric 
current. 


Fig'iVf  5.  Plots  of  coubtiint  c  luagiK-tir  lii-iil 


Figure  7.  Same  zis  figure  6  at  f  =  25.  top  left.  Maximum  magnetic  field  =  1.50. 
top  right.  Maximum  velocity  =  1.25.  bottom  left.  Maximum  contour  =  9.25. 
minimum  contour  =  -14.6  bottom  right.  Ma.ximum  contour  =  7.32.  minimum 


contour  =  -6.63. 


Figure  8.  Seime  as  figure  6  at  t  =  53.  top  left.  Maximum  magnetic  field  =  1.03. 
top  right.  Majcimum  velocity  =  0.36.  bottom  left.  Maximum  contour  =  1.S7. 
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minimum  contour  =  -1.S7  bottom  right.  Maximum  contour  =  1.92,  minimum 


X 

Figure  9.  Contour  plots  in  {x,y,z  =  0)  plane,  a.  Alignment  cosine  at  t  =  0.5. 
maximum  contour  =  0.910,  minimum  contour  =  -0.52S.  A  magnitude  of  1 
indicates  that  the  magnetic  field  is  aligned  with  the  electric  current  density. 
6.  a  at  t  =  0.5.  maximum  contour  =  2.51,  minimum  contour  =  -5.19.  c. 
Alignment  cosine  at  f  =  53.  maximum  contour  =  0.66S,  minimum  contour  = 
-0.879.  d.  o  at  <  =  53.  maximum  contour  =  4.62,  minimum  contour  =  -2.2S. 
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